function [output] = LOCA_main(input_array)

% data input
method = input_array.method;
data = input_array.data;
group = input_array.group_location;
alt_vec = input_array.alt_vec;

D = data.D;
X = data.X;
Y = data.Y;
        
switch lower(method)
    case 'ols'

        X_mat = [D X];
        betaHat_ols = X_mat\Y;
        resid = Y-X_mat*betaHat_ols;
        iXX = inv(X_mat'*X_mat);
        beta0_hat = betaHat_ols(1); % OLS estimate

        sctmp = cluster_se(X_mat,resid,iXX,group);
        se_cl = sctmp(1); % clustered standard error
        
        alt_power_curve = linspace(alt_vec(1),alt_vec(end),50);
        
    case 'iv'
        
        Z = data.Z;
        
        X_mat = [D X];
        Z_mat = [Z X];

        beta_iv = (Z_mat'*X_mat)\(Z_mat'*Y);
        beta0_hat = beta_iv(1);
        resid = Y-X_mat*beta_iv;
        iZX = inv(Z_mat'*X_mat);
        sctmp = cluster_se(Z_mat,resid,iZX,group);
        se_cl = sctmp(1); 
        
        alt_power_curve = linspace(alt_vec(2),alt_vec(end-1),50);
end

% p-value under the null
t_stat = (beta0_hat-0)/se_cl;
p_null = 2*(1-tcdf(abs(t_stat),max(group)-1));

% p-value under the alternatives
t_stat_alt = (beta0_hat+alt_vec-0)/se_cl;
p_alt = 2*(1-tcdf(abs(t_stat_alt),max(group)-1));

% p-value used for power curve
% alt_power_curve = linspace(alt_vec(1),alt_vec(end),50);
t_stat_power_curve = (beta0_hat+alt_power_curve-0)/se_cl;
p_power_curve = 2*(1-tcdf(abs(t_stat_power_curve),max(group)-1));


% output 
output.estimate = beta0_hat;
output.p_null = p_null;
output.p_alt = p_alt;
output.se = se_cl;
output.alt_power_curve = alt_power_curve;
output.p_power_curve = p_power_curve;



